##########################################################################################

library(data.table)
library(optparse)
library(dplyr)
library(ggplot2)
library(ggsci)
library(RColorBrewer)
library(ggpubr)

##########################################################################################

option_list <- list(
    make_option(c("--njmu_file"), type = "character") ,
    make_option(c("--tcga_file"), type = "character") ,
    make_option(c("--oncoSG_file"), type = "character") ,
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    njmu_file <- "~/20220915_gastric_multiple/dna_combine/baseTable/STAD_Info.addBurden.tsv"
    tcga_file <- "~/20220915_gastric_multiple/dna_combine/public_ref/TCGA/TCGA_STAD.TMB.tsv"
    oncoSG_file <- "~/20220915_gastric_multiple/dna_combine/public_ref/OncoSG/OncoSG_STAD.TMB.tsv"
    out_path <- "~/20220915_gastric_multiple/dna_combine/images/CompareBurden"
}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

njmu_file <- opt$njmu_file
tcga_file <- opt$tcga_file
oncoSG_file <- opt$oncoSG_file
out_path <- opt$out_path

dir.create(out_path , recursive = T)

###########################################################################################

col <- c(
  brewer.pal(9,"YlGnBu")[6],
  rgb(234,106,79,alpha=255,maxColorValue=255),
  rgb(203,24,30,alpha=255,maxColorValue=255),
  rgb(255,0,0,alpha=255,maxColorValue=255)
  )

names(col) <- c("IM" , "IGC" , "DGC" , "GC")

Variant_Types <- c("Missense_Mutation","Nonsense_Mutation","Frame_Shift_Ins","Frame_Shift_Del","In_Frame_Ins","In_Frame_Del","Splice_Site","Nonstop_Mutation")

###########################################################################################

dat_njmu <- data.frame(fread(njmu_file , quote=""))
dat_tcga <- data.frame(fread(tcga_file , quote=""))
dat_oncosg <- data.frame(fread(oncoSG_file , quote=""))

###########################################################################################
## 去除同一个患者多个样本 
dat_njmu <- subset( dat_njmu , Type!="IM + IGC + DGC")
dat_njmu <- dat_njmu[,c("Patient" , "Class" , "BurdenExon" )]
dat_njmu$From <- "NJMU"
dat_njmu <- subset( dat_njmu , Class %in% c("IGC" , "DGC") )
## 同一个人的合并
dat_njmu <- dat_njmu %>%
group_by( Patient , Class , From ) %>%
summarize( BurdenExon = median(BurdenExon) )
colnames(dat_njmu) <- c("Tumor" , "Class" , "From" , "BurdenExon")

###########################################################################################
## TCGA数据整理
dat_tcga$Class <- "DGC"
dat_tcga$Class[grep( "Intestinal" , dat_tcga$histological_type)] <- "IGC"
dat_tcga$BurdenExon <- dat_tcga$TMB..nonsynonymous.
dat_tcga$From <- "TCGA"
dat_tcga <- dat_tcga[,c("bcr_patient_barcode" , "Class" , "From" , "BurdenExon")]
colnames(dat_tcga) <- c("Tumor" , "Class" , "From" , "BurdenExon")

###########################################################################################
## Oncosg
dat_oncosg$Class <- "DGC"
dat_oncosg$Class[grep( "Intestinal" , dat_oncosg$Laurens.classification)] <- "IGC"
dat_oncosg$BurdenExon <- dat_oncosg$TMB..nonsynonymous.
dat_oncosg$From <- "EastAsian"
dat_oncosg <- dat_oncosg[,c("Patient.ID" , "Class" , "From" , "BurdenExon")]
colnames(dat_oncosg) <- c("Tumor" , "Class" , "From" , "BurdenExon")

###########################################################################################

dat_combine <- rbind( dat_njmu , dat_oncosg , dat_tcga )

###########################################################################################

my_comparisons <- list( c("EastAsian", "NJMU"), c("EastAsian", "TCGA"), c("NJMU", "TCGA") )
dat_combine$From <- factor( dat_combine$From  , levels = c("NJMU" , "EastAsian" , "TCGA" ) , order = T )

out_name <- paste0(out_path , "/MutationBurden.compare.tsv")
write.table( dat_combine , out_name , row.names = F , quote = F , sep = "\t" )

outcompare <- dat_combine %>%
group_by( Class , From ) %>%
summarize( BurdenExon = median(BurdenExon) )
out_name <- paste0(out_path , "/MutationBurden.compare.plot.tsv")
write.table( outcompare , out_name , row.names = F , quote = F , sep = "\t" )

###########################################################################################


plot <- ggplot( data = dat_combine , mapping = aes( x = From , y = BurdenExon , fill = From ))+
    geom_boxplot(alpha =1 , size = 0.9 , width = 0.6 , outlier.shape = NA) +
    facet_grid(.~Class) + 
    #geom_jitter(position=position_jitter(0.2) , aes(colour=From)) +
    scale_fill_npg() +
    scale_color_npg() +
    xlab(NULL) +
    ylab('Mutation rate per MB')+
    theme_bw() +
    stat_compare_means(comparisons = my_comparisons , label = "p.format") +
    #stat_compare_means(label.y = 1.1) +
    theme(panel.background = element_blank(),#设置背影为白色#清除网格线
        legend.position ='right',
        legend.title = element_blank() ,
        panel.grid.major=element_line(colour=NA),
        legend.text = element_text(size = 8,color="black",face='bold'),
        axis.text.x = element_text(size = 10,color="black",face='bold'),
        axis.text.y = element_text(size = 15,color="black",face='bold'),
        axis.title.x = element_text(size = 18,color="black",face='bold'),
        axis.title.y = element_text(size = 18,color="black",face='bold'),
        axis.line = element_line(size = 0.5)) 

out_name <- paste0(out_path , "/MutationBurden.compare.From.pdf")
ggsave(file=out_name,plot=plot,width=6,height=6)

###########################################################################################
my_comparisons <- list( c("IGC", "DGC"))

plot <- ggplot( data = dat_combine , mapping = aes( x = Class , y = BurdenExon , color = Class ))+
        geom_boxplot(alpha =1 , size = 0.9 , width = 0.6 , outlier.shape = NA) +
        geom_jitter(position = position_jitter(0.17) , size = 1 , alpha = 0.7) +
        facet_grid(.~From) +
        scale_color_manual(values=col) +
        xlab(NULL) +
        #ylim(0,15) +
        #scale_y_continuous(trans="sqrt" , breaks = c(0,0.5,1,2,3,4,5,10,20,30,40,50) ) +
        ylab("Mutation rate per MB")+
        theme_bw() +
        stat_compare_means(comparisons = my_comparisons,method = "wilcox.test") +
        theme(panel.background = element_blank(),#设置背影为白色#清除网格线
            legend.position ='left',
            legend.title = element_blank() ,
            panel.grid.major=element_line(colour=NA),
            plot.title = element_text(size = 8,color="black",face='bold'),
            legend.text = element_text(size = 10,color="black",face='bold'),
            axis.text.y = element_text(size = 10,color="black",face='bold'),
            axis.title.x = element_text(size = 10,color="black",face='bold'),
            axis.title.y = element_text(size = 10,color="black",face='bold'),
            axis.ticks.x = element_blank(),
            axis.text.x = element_text(size = 10,color="black",face='bold') ,
            axis.line = element_line(size = 0.5)) 
    
out_name <- paste0(out_path , "/MutationBurden.compare.IGC_DGC.pdf")
ggsave(file=out_name,plot=plot,width=8,height=6)

###########################################################################################
my_comparisons <- list( c("IGC", "DGC"))

dat_plot <- subset( dat_combine , From != "NJMU" )
dat_plot$Class <- factor( dat_plot$Class , levels = c("IGC" , "DGC") , order = T )

plot <- ggplot( data = dat_plot , mapping = aes( x = Class , y = BurdenExon , color = Class ))+
        geom_boxplot(alpha =1 , size = 0.9 , width = 0.6 , outlier.shape = NA) +
        geom_jitter(position = position_jitter(0.17) , size = 1 , alpha = 0.7) +
        facet_grid(.~From) +
        scale_color_manual(values=col) +
        #ylim(0,15) +
        xlab(NULL) +
        #scale_y_continuous(trans="sqrt" , breaks = c(0,0.5,1,2,3,4,5,10,20,30,40,50) ) +
        ylab("Mutation rate per MB")+
        theme_bw() +
        stat_compare_means(comparisons = my_comparisons,method = "wilcox.test") +
        theme(panel.background = element_blank(),#设置背影为白色#清除网格线
            legend.position ='left',
            legend.title = element_blank() ,
            panel.grid.major=element_line(colour=NA),
            plot.title = element_text(size = 8,color="black",face='bold'),
            legend.text = element_text(size = 10,color="black",face='bold'),
            axis.text.y = element_text(size = 10,color="black",face='bold'),
            axis.title.x = element_text(size = 10,color="black",face='bold'),
            axis.title.y = element_text(size = 10,color="black",face='bold'),
            axis.ticks.x = element_blank(),
            axis.text.x = element_text(size = 10,color="black",face='bold') ,
            axis.line = element_line(size = 0.5)) 
    
out_name <- paste0(out_path , "/MutationBurden.compare.IGC_DGC.noNJMU.pdf")
ggsave(file=out_name,plot=plot,width=7,height=6)
